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I. INTRODUCTION 


The dynamics of colloid suspensions is crucially influenced by flow-mediated correla¬ 
tions 1 2 . While these hydrodynamic interactions (HI) have an important role in the dynamics 
of ambient suspensions at thermal equilibriumP, their effect becomes even more pronounced 
for objects driven out of equilibrium, where the total force acting on each object generates a 
long-ranged flow, decaying as 1 /R with the distance R between the objects. A well-known 
example is colloid sedimentation, where HI lead to strongly correlated motions and large- 
scale dynamic structures 3 . Various types of driving, such as electrophoresis, are widely used 
to control the transport of colloids and other polyatomic objects 2 . Theoretical studies of 
driven colloids traditionally focus on regular particle shapes such as uniform spheres and el¬ 
lipsoids. The driving of more asymmetric objects is richeiEKS as it generally includes coupling 
between translation and rotation — when the object is subjected to a force it also rotates, 
and when it is under torque it also translated. The choice of a rotation sense under a 
unidirectional force implies a chiral response of the driven object. Such richer responses can 
be exploited to obtain “steerable colloids” — objects whose orientation and transport can 
be controlled in much more detail. For example, applying a torque by a rotating uniform 
magnetic field was used to achieve efficient transport of chiral magnetic objects 8 . Another 
example, which is the main issue of the present work, is the ability to achieve orientational 
alignment of asymmetric objects by applying an external forcd®d 

The earlier theoretical works of Refs. 8 11 dealt with isolated asymmetric objects in Stokes 
flow, which exhibit a chiral response. The object’s chiral response is encoded in the off- 
diagonal block of its self-mobility matrix, referred to as the twist matrix. Some objects 
have a twist matrix that leads them to align one axis in the body with the applied force. 
If the twist matrix has only a single real eigenvalue, the object becomes “axially aligned” 
in this way®^, and the aligning direction is along the corresponding eigenvector. Hence, in 
the absence of HI and thermal fluctuations, a set of identical, axially aligning objects reach 
a partially aligned state, where all the objects rotate about the same axis with the same 
angular velocity, but with an arbitrary phase. Furthermore, it was shown that, by applying 
an appropriate time-dependent forcing, the system can be driven to a fully aligned state, 
where all the objects are phase-locked with the force and rotate in synchrony 10 11 . 

In view of the above we use throughout this article the following terminology concerning 
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the response of various objects: (i) symmetric objects (such as a uniform sphere); (ii) regular 
objects, which are asymmetric objects with a vanishing twist matrix (such as a uniform 
ellipsoid); (iii) irregular objects, having a non-vanishing twist matrix; (iv) axially alignable 
objects, which are irregular objects, whose twist matrix has a single real eigenvalue. We 
note that the twist matrix depends on the position of the forcing point as well. For example, 
an ellipsoid whose forcing point is displaced from its centroid, i.e., an ellipsoid with a non- 
uniform mass distribution under gravity, has a non-vanishing twist matrix, and generally 
might be alignable. 

The theoretical groundwork for treating the HI between arbitrary objects in Stokes flow 
was laid by Brenner and O’Neill 1213 . The theory was subsequently applied to a pair of 
particles of various regular shapes 11 20 . To this one should add many earlier studies of 
the collective dynamics of suspensions made of ellipsoids 21 25 . We note that there are key 
differences between asymmetric objects, such as ellipsoids, and the irregular objects studied 
here. The symmetries of a uniform ellipsoid lead to: (a) the absence of a translation- 
rotation coupling for a single object, and therefore lack of alignability; (b) the absence of 
a 1 /R 2 contribution to the relative velocity developed between two such objects at mutual 
distance R. Finally, several numerical techniques have been introduced to treat suspensions 
of arbitrarily shaped object^®®. 

In this work we focus on simple, general properties of the pair HI between two arbitrarily 
shaped objects at zero Reynolds number, and the resulting effect on their orientational 
alignment. The study of translational effects will be presented in a separate publication. 

The work is made of two distinct parts. The first part treats rigorously the instanta¬ 
neous hydrodynamic interaction, i.e., the pair-mobility matrix. We use Brenner’s analytical 
framework 31 22, specializing to the leading order of the HI in the distance between the objects 
(multipole expansion, also known as the method of reflection JD). The second part addresses 
the time-dependent trajectories of forced objects. This is a multi-variable, highly non-linear 
dynamical system exhibiting complex and diverse dynamics. In this part we are limited to 
numerical integration of the objects’ trajectories. We provide typical examples for the time 
evolution of pairs of stokeslet objects. 

We begin by discussing in Sec. [TT] the general properties and symmetries of the pair- 


mobility matrix for two arbitrarily shaped objects. In Sec. HI we apply a multipole expansion 
to the pair-mobility matrix and obtain results for the instantaneous HI at large distances. 
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In Sec. [TVjwe derive the resulting properties of stokeslet objects, and in Sec. [V]we use them 
to perform numerical time integration for the evolution of object pairs and their alignment. 


Finally, in Sec. VI. we discuss several consequences of our results. 


II. PAIR-MOBILITY MATRIX: GENERAL CONSIDERATIONS 
A. Structure of the Pair-Mobility Matrix 


The kinematics of a rigid object is represented by a translational velocity V, which refers 
to an arbitrary reference point rigidly affixed to the object, and an angular velocity uj. We 
designate the reference point as the origin of the object. Note that the angular velocity of 
the object is independent of the choice of its origin, and that the origin does not necessarily 
lie on the instantaneous axis of rotation of the object. 

Consider two arbitrarily shaped rigid objects, a and b, with typical size l, subject to 
external forces and torques F a , F b and r a , in an unbounded, otherwise quiescent fluid 
of viscosity rj. In the creeping flow regime, the objects respond with linear and angular 
velocities to the external forces and torques through a 12 x 12 pair-mobility matrix , 



( 1 ) 


where we define generalized velocity and generalized force 6-vectors, V x = (V x , Iuj x ) t and 
jz x — (F x ,t x /1) t for x = a, b. The diagonal blocks, M aa and M bb , correspond to the self¬ 
mobilities of the objects (which nevertheless depend on the configuration of both objects). 
The off-diagonal blocks, M ab and M ba , describe the pair hydrodynamic interaction. We 
hereafter omit the factor (■ rjl ) -1 (i.e., set rjl = 1). This, together with the representation of 
the generalized forces and velocities, make M dimensionless and dependent on the geometry 
alone. Throughout the text we designate 6-vectors and matrices with calligraphic font and 
blackboard-bold letters, respectively. A detailed description of the notation used in the 
article is given in Appendix [A] 

Since V and t depend on the choice of object origins, so does the pair-mobility matrix. 
The transformation between pair-mobility matrices corresponding to different origins is given 
in Appendix |Bj 
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The pair-mobility matrix is a function of the objects’ geometries, their orientations, and 
the vector connecting their origins, indicated hereafter by R. (We define the direction of R 
from the origin of object b to the origin of object a.) The geometry of object x is denoted 


by r x . For example, if the object consists of a discrete set of N x stokeslets (see Sec. IV A), 
then r x is a 3iV x -vector specifying the positions of the stokeslets; otherwise, it represents 
the surface of the object. 

The pair-mobility matrix is positive-definite and symmetric 13331 . Hence, M a6 = (M fea ) T , 
and the self-blocks can be written as 






As in the analysis for isolated object^®, the self-mobility matrix contains the following 3x3 
blocks: the alacrity matrix A (translational response to force); the screw matrix § (rotational 
response to torque); and the twist matrix T (translation-rotation coupling). The twist 
matrix characterizes the chiral response of the object (the sense of rotation under a force). 
In the present article we deal with alignable objects, whose individual T is necessarily non¬ 
vanishing. Furthermore, in the case of a pair of objects, the presence of the other object 
makes the self-twist matrix, T xx , differ from the single-object one. As to the off-diagonal 
blocks of the pair-mobility matrix, the symmetry of M implies the following structure: 


M ab = I 

f A ab 

(' J ba ) T ) 

1 , M ba = I 

f (A ab ) T 

(J ab ) T \ 

! 

[ TP 6 

S ab j 

1 1 

i T 6a 

(§ a6 ) T ) 


B. Further Symmetries of the Pair-Mobility Matrix 

The discussion in the preceding subsection has been for a general pair of objects, which 
are not necessarily identical. In the present subsection, we focus on the case in which the two 
objects are identical in shape and orientation , i.e., r a = r b = r. Our goal is to understand 
what the instantaneous relative velocities (linear and angular) between the two objects are, 
when the objects are subjected to the same external forcing. The restriction to identical 
objects makes M invariant under exchange of objects. This additional symmetry is made of 
two operations: interchanging the blocks M aa M bb and M afe M fea ; and inversion of R. 
That is, 

M(r, R) = EM(r, -R) E" 1 , (2) 
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where E is a 12 x 12 matrix which interchanges the objects, 


E = 




with Igx6 denoting the 6x6 identity matrix. 

The symmetry to object exchange, when combined with the parity of M (i.e., whether 
it remains the same or changes sign) under i?-inversion]— has important consequences for 
the effect of hydrodynamic interactions on alignment. If M has a definite parity one can 
determine what the relative response of the objects to forcing is — i.e., whether they attain 
the same or the opposite linear and angular velocities. If the term is symmetric to inversion, 
the velocities would be identical, and if it is antisymmetric, they would be opposite. This is 
because 


/ M aa (R) M a6 (i?)\ f M aa (-R) M afe (—i?)\ /M bb (R) M ba (R) 

l M ba (R) M bb (R) J l M ba (-R) M bb (—R) J l M ab (R) M aa (R) 

where the second equality comes from the response to exchange of objects, Eq. ([2]). Conse¬ 
quently, under identical forcing of the two objects one finds, 



V a = (M aa + M ab ) f = ± (M“ + M ba ) f = ±V b . 


( 4 ) 


Thus, since any M can be decomposed into even and odd terms, we find that only the odd 
ones cause relative motions of the two objects. 

The pair-mobility as a whole, however, never has a definite parity under i?-inversion, 
i.e., it is made of both even and odd terms. This becomes clear when M(r, R) is expanded 
in small l/R, i.e., in multipoles. A general discussion of the parity of each multipole term 
is given in the next section. For now, let us consider those two leading multipoles which 
are independent of the objects’ shape, and therefore always exist. The monopole-monopole 
interaction (Oseen tensor), which is the leading term in A ab making particle a translate due 
to the force on particle b, is symmetric under /^-inversion. The part of the monopole-dipole 
interaction causing the second object to rotate due to the force on the first, i.e., the leading 
term in T a6 , is antisymmetric. For example, even the most symmetric pair of objects — two 
spheres — has an A-symmetric A ab , leading to zero relative velocity, and an A-antisymmetric 
T afc , causing them to rotate with opposite senses 1 . Thus, for a general object, the highest 
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order which maintains M of definite parity is the monopole 1/R Oseen one, which is even. 
(The self-blocks are constant up to order 1/R 4 ] see below.) 

From this discussion we can immediately conclude that, to leading order in the separation 
of two identical, fully aligned objects, their instantaneous hydrodynamic interaction must 
linearly degrade the alignment. The leading degrading term comes from T ab , their rotational 
response to force, and is of order 1/R 2 . ft is worthwhile to note again that such a rotational 
response is present as well for a pair of uniform spheres or ellipsoids; yet, such regular objects 
are not alignable to begin with. 

The relation between object-exchange symmetry and the symmetry of the linear-velocity 
response is intimately related to the issue of hydrodynamic pseudo-potential^, which will 
be discussed in detail in a forthcoming publication. 

III. FAR-FIELD INTERACTION: MULTIPOLE EXPANSION 

There are two characteristic length scales in our problem: the typical size of the objects, 
l, and the distance between them, R = |f?|. If l <C R, we can write the pair-mobility matrix 
as a power series in ( l/R ), 


M = M ( o) + M ( i) + M(2) + ..., 


where M( n ) ~ ( l/R ) n . The analysis of this expansion as given below holds for any pair 
of objects, whether identical or not. The zeroth order, M(o), is a block diagonal matrix 
which is made of the self-mobilities of the two non-interacting objects. (These should be 
distinguished from M aa and M 66 , the self-mobilities of the interacting objects.) 

The hydrodynamic multipole expansion (also known as the method of reflections) is based 
on the Green’s function of Stokes flow, the Oseen tensor 1 } given in our units (■rjl = 1) by 



1 l 
8tt r 



(5) 


which is a symmetric 3x3 tensor, invariant under F-i rivers ion. A point force at Fo, S(f—fo)f, 
generates a velocity field u(f) = G(F— Fo) • /. 

We obtain two general results concerning the multipoles of the hydrodynamic interaction 
between two arbitrary objects. The two objects need not be identical. The proofs are given 
in Appendix [Dj 3 - 
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1. The leading interaction multipole in the self-blocks of the pair-mobility matrix is n = 4. 
That is, any response of one object to forces on itself, owing to the other object, must 
fall off with distance R between the objects at least as fast as R~ 4 . 

2. The nth multipole has self-blocks of (— l) n parity, and coupling blocks of the opposite, 
(_l)n+i p ar ity. Thus, e.g., the leading term in M aa , proportional to R~ 4 , is invariant 
under R-inversion, and the R~ 4 part o/M afe changes sign under R-inversion. Likewise 
for the multipole varying as R~ 5 , the M aa changes sign under R-inversion while M a6 
remains invariant , 35 


These statements pertain to the mobility matrix. As to the propulsion matrix (the inverse 
of the mobility matrix), the leading correction to the self-block becomes ~ 1 /R 2 , and the 
second statement concerning parity remains intact. 

We now consider for a moment two identical objects and specialize to the first and 
second multipoles, i.e., the hydrodynamic interaction up to order 1/R 2 . The discussion in 
the preceding and current sections implies the following form of the two leading terms in 
the pair-mobility matrix: 



( 6 ) 


In more detail: there are no first- and second-order corrections to the objects’ self-mobility. 
Hence, these two multipoles have definite parities — the first is even, and the second is odd. 
Consequently, the first multipole does not cause any relative motion of the two objects, 
whereas the second mutipole makes them translate and rotate in opposite linear and angular 
velocities. 

The essential characteristics of the first two multipoles are schematically illustrated in 
Fig. B The first multipole arises directly from the Green’s function, 




(7) 


where G (R) is the Oseen tensor, given in Eq. ((sj) . 

In the interaction described by the second multipole one object sees the other as a point, 
see Fig. [Tj Accordingly, this term contains two types of interaction: (1) the response of 
object a to the non-uniformity of the flow due to the force monopole at object b (regarded 


as a point); (2) the advection of object a (regarded as a point) by the flow due to the 
force dipole acting at object b. These two effects are both proportional to VG(-R) ~ 1/R 2 . 
Each can be written as a product of a tensor which arises from the medium alone, through 
derivatives of the Oseen tensor VG(R), and another tensor which depends on the objects’ 
geometry. The second-order correction to the velocity of object a is given by the sum of 
these two effects, each expressed in terms of a coupling tensor 0 and an object tensor <f> 


V ( °2) = Mg • P 

Mg = $ a :Q(R)~ O t (R) : P, (8) 


where the double dot notation denotes a contraction over two indices. Equation (|8]) contains 
three tensors of rank 3, denoted by capital Greek letters. The first, $, with dimensions 
6x3x3, gives the generalized velocity of the object in linear response to the velocity 
gradient of the flow in which it is embedded. The second, <f>, having dimensions 3x3x6, 
gives the force dipole acting on the fluid around the object’s origin in linear response to the 
generalized force acting on it. Both <f> and <f> depend on the objects’ geometry alone 38 . The 
third tensor, 0, with dimensions 3x3x6, describes the coupling of these object responses 
through the fluid. It is given by 


Qskj(R) 


d s G kj {f)\^ j = 1,2,3 
0 j = 4,5,6. 


(9) 


Repeating the same procedure for V b in response to T a while using the odd parity of 0, we 
get 

Mg = G t (R) : P - P : G(R). (10) 


The tensors and $ are not independent 39 . We now show that $ = < h T . The symmetry 
of M implies that each multipole is also a symmetric matrix. Using Eqs. (|8|) and (10) and 
equating (Mg) 7 ’ = Mg, we get <3> a = (<f> a ) r and <f> b = ($ fe ) T 
To summarize, the matrix M/ 2 ) is given by 



( P _ 

\-<f> h : e(R) + [$ a : 0(R)] T 


: Q(R) - [$ b : 0(R)] T 
0 


( 11 ) 


This results is valid for a general pair of objects. If the two objects are identical, the off- 
diagonal blocks have the same form with opposite signs. The additional condition that the 
entire M must be symmetrical implies then that each block by itself is antisymmetric. 
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Mg } - VG (R) ~ 


FIG. 1. Illustration of the two leading orders of the hydrodynamic interaction between two forced 
objects. The leading term in the pair-mobility matrix (light blue/dash-dotted arrow between the 
objects’ origins), decaying as 1/R, comes from the point-like response of object a to the local flow 
caused by the force monopole on object b (blue/thick arrow). The next-order term, decaying as 
l/R 2 , has two contributions: (i) The point-like response of object a to the local flow caused by the 
force dipole on object b (red/dashed arrow from the red/thin arrows at object b to the origin of 
a), (ii) The response of object a to the local flow gradient caused by the force monopole on object 
b (magenta/dotted arrow from the origin of b to the magenta/thin arrows at object a). 


By separating the tensors <3> and 0 into their symmetric and antisymmetric parts, the 
second-order term of the pair-mobility matrix can be simplified further. It should be men¬ 
tioned, in addition, that the d> tensor depends on the origin selected for the object. These 
two technical issues are addressed in Appendices [E] and [Cj respectively. Finally, we note that 
the terms in these tensors corresponding to the translational response vanish for spheres and 
ellipsoids. Consequently, two such regular objects develop relative velocity only to orders 
1/R 3 and above. 
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IV. NUMERICAL ANALYSIS FOR STOKESLET OBJECTS 


In the preceding sections we have derived the general properties of the instantaneous 
hydrodynamic interaction between two arbitrarily shaped objects. We now move on to the 
second part of the work, addressing the time evolution of the two objects. This complicated 
problem is not tractable analytically, and we resort to numerical integration of specific 
examples. Because of the complexity of the problem, and since we are interested in generic 
properties, we allow ourselves to restrict the analysis to the simplest, even if unrealistic, 
objects. Arguably the simplest form of an arbitrarily shaped object is the so-called stokeslet 
object —a discrete set of small spheres, separated by much larger, rigid distances, where 
each sphere is approximated as a point force. The sparseness of these objects makes them 
free-draining, which may be valid for macromolecules but not for compact objects. 

We treat pairs of identical objects, each made of four stokeslets. To obtain representative 
sampling of numerical examples we do not design these objects but create them randomly. 
Four points are placed at random distances ranging between 0 and 1 from an arbitrary origin. 
The origin is then shifted to the points’ center of mass. The radius p of the stokeslets is taken 
as 0.01. The resulting configuration is checked to be “sufficiently chiral”, in the sense that 
the T-matrix of the individual object is strongly asymmetric, having a single real eigenvalue 
of absolute value |A 3 | > 0.005, which makes the object axially alignable. (See Sec. [TJ). 
Examples of the stokeslet objects we use are provided in Fig. [2j 

The way to calculate the mobility of a single stokeslet object was presented in Ref. Sj 


First, we briefly present in Sec. IV A the simple extension of this method to pair-mobilities. 
We calculate both the pair mobility and the tensor <F introduced in Secs.[TT]and III The latter 
allows us to calculate pair mobilities up to second order in the multipole expansion. Section 


IV B describes how we use the pair mobility to numerically calculate the time evolution of 


the pair configuration. 


A. Pair-Mobility and $ Tensor 

The properties of a stokeslet object can be derived self-consistently from the linear rela¬ 
tions which describe the stokeslets’ configuration. This is done without finding the stokeslets’ 
strengths explicitly. Below we fold the pair-mobility matrix, and the tensor associated 
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with a single object, given the stokeslet configuration and the size of the spheres that they 
represent. 

Each of the two objects, x = a, b, consists of N x stokeslets, F x = (Fj x ,..., ) , in a 

known configuration, r x = Here, we use the notation of a bold letter to 

denote a set of N 3-vectors, and r x indicates the position 3-vector of the nth stokeslet in 
object x with respect to the object’s origin. Each stokeslet is a sphere of radius p, where 
p < min(rf,..., r x N ). The boundary conditions at the sphere surface enter only through its 


self-mobility coefficient. The velocities of the spheres, are known from the object’s linear 
and angular velocities, 




^3X3, -rT x /A 


with U x = 


for x = a,b, 


( 12 ) 


w ' ' ' 7 \W-f5„7 1 ) 

where the matrix y x obtained from the vector y is defined as ( y x )ij = e^yk- Each stokeslet 
force is proportional to the relative velocity of the sphere that it represents, with respect to 
the flow around it as created by the other stokeslets. This gives a linear relation between 
the stokeslets and the velocities of the spheres 4 ^, 



(13) 


(L: 


XX \ 
nm.) ij 


^*3 

7 -1 £ 


(Tn - r£) if n ± m 




else 


y(i2 + r^-f*), (14) 

and 7 = 67T p/l. 

First we find the pair-mobility matrix as a generalization of the analysis in Ref. [9] The 
sum of the stokeslets and the corresponding total torque must be equal to the external 
generalized forces applied on the objects. In a matrix form we can write 

(15) 


V 

l-( 


! \ 



Using Eqs. (12) and (13), we have 

T 

L aa h ab ' 
h abT L bb 
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From this expression we identify the pair-mobility matrix as 



L aa h ab 
h abT h bb 



-l 


( 16 ) 


This expression allows to calculate the pair-mobility matrix, with the help of Eqs. (12) and 


(14), based on the stokeslets’ configuration and the Oseen tensor alone. 


Next, we derive the <h‘ T tensor of a stokeslet object x. From this tensor we may readily 


obtain the second multipole of the pair interaction (cf. Sec. III). The force dipole around 


the origin of a forced object is given by [Eq. (|8j)] , ( rF) x = (<3R) T • T x . Similar to the U x 
matrix relating the stokeslets to the total generalized force, F x = (U x ) r • F x , we define 
a tensor of rank 3, T x , which relates the stokeslet forces to the total force dipole on the 
object by ( rF) x = ( T x ) t ■ F x . (Note that no force dipole is applied on the individual 
stokeslets; being arbitrarily small they possess only a force monopole.) Specifically, it is 
made of N blocks of 3 x 3 x 3, given by (T n )ij S = r ntS Sij, n = 1... N, i,j,s = 1, 2, 3 (i.e., 


r n:S is the s Cartesian coordinate of the stokeslet n). Using Eqs. (12) and (13), we have 
( rF) a = (T“) T ■ (Ld®)" 1 ■ U a ■ V a . This implies (<h*) r = ( T x ) t ■ (LT 1 ) -1 • IF ■ M* lf . Recalling 
that the matrices M se if and L are symmetric, we finally get 


cjp = 


p 

lelf 


(IF) T • (L^)" 1 ■ T x . 


(17) 


It is important to note that in the above derivation we compute M and under the 
assumption that, for each object, the stokeslet sizes are arbitrarily small compared to the 
distances between them, p -C l (where l is the object’s radius of gyration). However, in a 
more general case one can use the Rotne-Prager-Yamakawa tensor^^!, which corrects the 
Oseen tensor for force distributions with finite size 28 . 


B. Numerical Time Integration 

We present a numerical integration scheme for the dynamics of two stokeslet objects. We 
should first define the reference frames used in the scheme. Each rigid object is characterized 
by axes which are affixed and rotate with it. We define the object reference frame (ORF) 
such that its z axis coincides with the object’s alignment axis (the corresponding eigenvector 
of the T-matrix). The other two axes are selected arbitrarily. The z axis of the laboratory 
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frame is defined along the forcing direction. During the evolution we follow the translation 
and rotation of the ORF in the laboratory frame. 

We represent the orientation of an object by the Euler-Rodrigues 4-parameterd 4 ®, defined 
by (r, D) = (cos h sin §), where n and 9 are the axis and angle of rotation 11 . The following 
properties hold for this 4-parameter representation: (a) The norm of (T, D) in 4D-space is 
unity, T 2 + D 2 = 1. (b) A rotation matrix is given by Rodrigues’ rotation formula, 

=I 3x3 + 2m x + 2(D X ) 2 . (18) 


(c) The parameters are invariant under inversion, i.e., (T,h2) and (—T, —O) correspond to 
the same orientation, (d) Given the angular velocity of the object, the dynamics of its 
orientation simply reads 


0 —cu T 


, U U! 



(19) 


Since we choose the ORF such that the z-axis is the axis of alignment, the terminal 
orientation of an axially alignable object under a constant force along the z-axis is (r, D) = 
(cos(^i^), z sin( ii ^ 1 ^)), where ex is a constant phase which depends on the object’s initial 
orientation at time t = 0. 

The state of a pair of objects at time t is described by the position of the origins of the 
objects, R a (t) and R b (t), and their orientation parameters, (T a (t),Cl a (t)) and (T b (t),Ct b (t)). 
We time-integrate from the initial state, Rq = (0, 0, 0), Rq — Rq = Rq, (Tg, Dq) and (Tq, Dq), 
as follows. Given the positions of the stokeslets at time t, the pair-mobility matrix, M(t), 
is calculated as explained in Sec. IV A either exactly or using the multipole approximation. 
Then, the linear and angular velocities of the objects are given by (V a (t), V b (t)) 7 = M(t) ■ 
(J- a (t),J- b (t)) T . From them the origins and orientations of the objects at time t + dt are 
derived according to 


R x (t + dt) = R x (t ) + V x (t)dt 


T x (t + dt) 
Q x (t + dt) 


= exp 



-u 


oT 


X T 



( 20 ) 

( 21 ) 


for x = a, b. During the evolution we make sure that the small stokeslet spheres do not 
overlap, and that the pair mobility matrix remains positive-definite. In practice we never 
encountered such problems when using the exact pair mobility matrices; when it did happen 
in the case of the multipole approximation we stopped the integration. 
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We define a scalar order parameter which characterizes the degree of mutual alignment 
of the two objects, 


m = 


(T a ,fi a ) • (T b ,Q b ) 


= r a r b + n a -n l 


( 22 ) 


As required, the order parameter is invariant under 3-rotation. This can be verified by 
explicitly applying a 3-rotation to the laboratory frame, or alternatively, by the following 
argument. Since 3-rotation leaves the norm of the 4-parameter orientation unchanged (prop¬ 
erty (a) above), it is a unitary transformation in 4-space. Hence, the dot product is invariant. 
When the objects are aligned, (T a ,f2 a ) = ±(T 6 ,f2 6 ), and m = 1; otherwise 0 < m < 1. In 
the case of partial alignment, m = cos 2 (4r), where A a is the mutual phase clifferencd^. 

Another scalar property of the two-object system is the energy dissipation rate. At time 
f, the latter is given by V a {t) ■J ra {t) + V b (t) ■ F b (t). Since the pair-mobility matrix is positive 
definite the energy dissipation of the driven pair is positive at all times. 


V. NUMERICAL RESULTS: EFFECT ON ALIGNMENT 

We present in Figs. [3]-[8] several examples for the numerically integrated evolution of 
object pairs under various conditions. One can be immediately appreciate the diversity 
of possible trajectories. To make your way through this richness it is important to make 
two distinctions between types of trajectories. The first distinction is between constant 
forcing (as in sedimentation), which can make the objects only partially aligned without 
synchronizing their phases of rotation 5 ®, and a time-dependent forcing protocol, which is 
known to lock the phase of an individual object onto that of the forc^ 011 . The main 
issue examined below is how the presence of hydrodynamic interaction affects these two 
behaviors. The second distinction, therefore, is whether hydrodynamic interactions are 
included (dashed, dotted and dash-dotted/colored curves) or turned off (solid gray curves). 
In the absence of hydrodynamic interactions (or when they get weak as the objects move far 
apart), the time-dependent aligning force will make the objects fully synchronized, whereas 
under constant forcing the objects will generally become unaligned. 

The results are presented in a dimensionless form, using units such that rj = |a>o| = 
1 and p = 0.01. The distances between the stokeslets of each object are taken ran¬ 
domly between 0 and 1; hence, p <C l ~ 1. The time-dependent forcing protocol is F = 
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F 0 (—sin(wot) sin( 6 *), cos((not) sin( 6 *), — 008 ( 6 *)), where 6 = 0.17T, F 0 = — |A 3 1, cno = sign(A 3 ) 
and A 3 is the real eigenvalue of the single-particle twist matrix. We examine both the 
trajectories of the separation vector connecting the origins of the two objects, and the cor¬ 
responding evolution of the orientation order parameter. 

We begin with the case of a time-dependent forcing, Figs. [3] and [4} The first observation, 
most clearly demonstrated in Fig. 0b), is that hydrodynamic interaction degrades the align¬ 
ment of the two objects, as has been rigorously inferred based on symmetry considerations 
in Sec. II B| Another conclusion, supported by additional examples not shown here, is that 
most objects, which start sufficiently far apart, especially if they start fully aligned, tend to 
repel each other (Fig. 0. Even if they are not fully aligned, the growing distance and weak¬ 
ening interaction make them individually more aligned with the forcing, and therefore also 
mutually synchronized. Thus, the repulsion helps restore the alignment at long times. The 
increasing separation occurs in the xy plane, while along the z axis the separation decreases 
and saturates to a finite distance, dependent on initial conditions, see Fig. [3j 

The repulsion is accompanied by a decrease in dissipation rate (up to small oscillations), 
as demonstrated in Fig. [7} When the HI is turned off, the dissipation rate reaches a constant 
value as the two independent objects set into their ultimate aligned state (solid curves in 

Fig.0. 

The repulsive effect is observed for most examples of our randomly generated pairs of 
objects but is not a general law. For instance, when the objects start at a sufficiently small 
separation, some pairs remain “bound” in a limit cycle, oscillating about a certain mean 
separation and mean orientational alignment, as demonstrated by the green/dashed curves 
in Figs. |3] and [4| 

In Figs. [5] and [ 6 ] we examine the same properties under constant forcing. The two effects— 
degradation in the alignment of a pair which is initially fully aligned, and mutual repulsion— 
are observed here as well. Yet, in the absence of a time-dependent aligning force, as the 
two objects move apart, alignment is not restored. At long times, and for the common 
case of repulsion, we distinguish between two observed behaviors: a) The order parameter 
continues to change without saturating to a constant value (e.g., red/dash-dotted curve and 
cyan/dotted curve in Figs. 0a) and 0 b), respectively). This non-intuitive result can be 
explained as follows. The fact that the interaction becomes weak does not necessarily imply 
that the accumulation of phase difference stops. If the two distant objects are partially 
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aligned we have m ~ cos 2 (8uj z t/2), where Suj z is the difference between the objects’ angular 
velocity along the alignment axis. Hence, if the decay of 8uj z with time is slower then f -1 
then phase difference will continue to accumulate. This depends on the detailed dynamics of 
repulsion which will be addressed in publication II. b) The other option is that m converges 
to some value dependent on the initial state, with no particular chosen m (green/dash-dotted 
curve in Fig. |b](b) and cyan/dash-dotted curve in Fig. |b|c)), i.e., the two objects continue to 
rotate with a fixed relative orientation. In the examples that we checked there seems to be a 
tendency toward ultimate anti-alignment (m = 0). Therefore, we also checked the stability 
of anti-alignment in pairs which start from such a state. Fig. |6^d) examines the stability of 
this configuration for objects initially confined to the xy plane (perpendicular to the force). 
Whereas the aligned pair (blue/dotted curve) is unstable, the anti-aligned one (red/dashed 
curve) remains stable for the duration of integration. It may well be that this stability 
survives for a long but finite time, see e.g., dark red/dotted curve in Fig. |6](c). In addition, a 
separation of the pair along the £-axis destabilizes an anti-aligned pair as well (examples not 
shown). Finally, we note that even if the final phase difference were arbitrary and uniformly 
distributed, the value of m would be evenly distributed around 1/2 but non-uniformly, with 
larger weights on m = 0,1. (This follows from the definition of m, see Eq.([22|).) 

Figure [8] compares results obtained using the full pair-mobility matrix of the stokeslet 
objects with those obtained from the multipole (dipole) approximation. As expected, the 
two calculations agree for objects whose mutual distance increases with time, and disagree 
for objects whose trajectories reach close proximity. 

Further investigation (not shown here) of the orientational dynamics of the objects sug¬ 
gests a possible explanation for the characteristic repulsion between two chiral objects. 
In the absence of HI, each object rotates along F and translates on average along F. 
One contribution to the dipolar term of the HI comes from the effect on each object by 
the vorticity of the Oseen flow caused by the other object. This perturbative angular 
velocity is along an axis which is perpendicular to the separation vector and the force, 
Wj ow oc —R x F and a)| ow oc Rx F. The competition between this rotation and the aligning 
self-response of each individual object results in an inclination of the two objects relative 
to their non-interacting state. This inclination alters the average unperturbed linear ve¬ 
locity of the object by a small rotation about the R x F direction— counter-clockwise for 
object a and clockwise for object b. Hence, the two objects glide away from each other, 


17 


2 N 



FIG. 2. Two examples of axially alignable stokeslet objects, which were used in the simulations. 
The objects comprise four stokeslets connected by dragless rods. The origin of the objects is at 
point (0,0,0) and the aligning direction is —z. The object on the left corresponds to the dark 
red/dotted trajectories in the left panels of Figs. [5] and [6j and the one on the right corresponds to 
the purple/dashed trajectories in the right panels of Figs. [3] and [4| 


R 2 = 2(V a — V b ) • R oc ((—R x F) x F) ■ R — R{ 1 — {Ft, ■ F ) 2 ) > 0, where the proportionality 
constant is positive, i.e., the separation increases with time (unless R || F. for which the 
whole argument does not hold). 


VI. DISCUSSION 

Irregular objects display rich dynamics already at the level of a pair of objects, as has been 
demonstrated above. In the present work we have focused on the effect of the hydrodynamic 
interaction on the orientational alignment of such object pairs. 

The hydrodynamic interaction, in general, degrades the alignment. We have rigorously 
proven the instantaneous linear degradation for fully aligned objects at large mutual dis¬ 
tances. In other circumstances, such as nearby or unaligned objects, the hydrodynamic 
interaction may have an opposite effect. The leading degradation effect in distance is dipo¬ 
lar rather than monopolar; yet, it is significant — a large mutual distance (compared to 
the object size) is required to make the degradation negligible. More quantitatively, the 
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FIG. 3. Trajectories of object separation under time-dependent forcing. The three rows, from top 
to bottom, correspond, respectively, to the separation along the z direction, its projection onto the 
xy plane, and its total magnitude. The squares in the middle row indicate the state at the end of 
the simulation. The panels show results for three different objects, starting from either a random 
mutual orientation (left column) or their fully aligned state (right column). The green/dashed 
trajectory on the right panels was integrated longer than 150 time units to verify that it continues 
in a limit cycle. 
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FIG. 4. Orientation order parameter as a function of time, under time-dependent forcing, for the 
examples of Fig. [3j (a) results for random initial orientations (examples on the left column of Fig. 
[3]); the additional solid gray curves correspond to non-interacting objects, (b) results for initially 
fully aligned object pairs (right column in Fig. [ 3 ]). 

degradation will be significant when the perturbation to the angular velocity due to HI, 
Soj, becomes comparable to the inverse of the time required to align a single object. The 
unperturbed angular velocity is given by loq = T se ifF. The dimensionless eigenvalue of the 
self-twist matrix is generally found to be about an order of magnitude smaller than the di¬ 
mensionless self-mobility coefficient 810 -^, i.e., uq ~ 0.1F/(87t/ 2 ). As presented in Sec. |II B[ 
5oj ~ T ab F ~ F/(87r/ 2 )(Z/i?) 2 . The alignment time is typically f a 1 ~ 10/ojq (see Fig [ 4 ]) . 
Hence, the degree of degradation is t a \Suj ~ 10 2 {l/R) 2 . The conclusion is that the separation 
between the objects should be larger than ten times their size to maintain alignment. In the 
case of many objects this implies a maximum volume fraction ( l/R ) 3 ~ 10 -3 . 

At the same time, as shown in Sec. [VJ for most of our randomly generated pairs of 
objects, the hydrodynamic interaction makes the rotating objects repel each other. As a 
result, at long times the hydrodynamic interaction usually becomes negligible and each of the 
objects gets aligned again with the time-dependent force. In that section we also provided 
a possible explanation for this repulsion, related to the mutual rotation of the two objects 
which causes them to glide away from each other. I 11 fact, the objects need not be irregular 
to exhibit this gliding effect; two forced ellipsoids which start parallel to one another will 
experience the same repulsion! 20 * 29 * 46 ! The resulting hydrodynamic “pseudo-potential” 364 ' 
will be addressed in a future publication. For the case of a finite number of objects the 
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FIG. 5. Trajectories of particle separation under constant forcing. The meaning of the various 
panels is the same as in Fig. |3j In all the examples shown here, the two objects repel each other 
except of the example which corresponds to the blue/dashed curve in the left panels. (The repulsive 
trajectories were actually integrated to times longer than presented here.). 
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FIG. 6. Orientation order parameter as a function of time, under constant forcing, for the examples 
of Fig. [5j (a) results for random initial orientations (examples on the left column of Fig. [5]); the 
solid gray curves correspond to non-interacting objects, (b) results for initially fully aligned object 
pairs (right column in Fig. |3j) . (c) results for objects with initial partial alignment (rotating around 
the same axis with random initial phases), (d) the stability of anti-alignment; shows trajectories 
of two identical pairs, which start on the xy plane from the same separation and axes of rotation 
but with different relative phases. Blue/dotted and red/dashed curves represent, respectively, a 
pair which starts aligned (zero relative phase) and one which starts anti-aligned (relative phase of 
Vr). 

repulsion will help restore the alignment as the objects drift apart. It should be kept in 
mind, however, that the repulsion is not a general law. We observed it for a few dozens 
pairs of stokeslet objects. As mentioned above, it also holds for a pair of well separated 
ellipsoids. Yet, a few counter-examples have been also provided in SecjV} 
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FIG. 7. Dissipation rate as a function of time for object pairs starting from arbitrary orientations, 
under time-dependent forcing (a) and constant forcing (b). Dash-dotted and dotted colored curves 
correspond to the examples of the same styles/colors in the preceding figures. Solid curves show 
the results in the absence of HI. 




FIG. 8. Comparison between the evolution of pair separations obtained using the full pair-mobility 
matrix (dashed, dotted and dash-dotted colored curves) and its multipole approximation (solid 
curves). Each panel presents three examples of pairs under time-dependent (a) and constant 
forcing (b). All pairs start from a fully aligned state. The multipole approximation includes the 
monopolar and dipolar terms. 

An interesting counterpart of the effects discussed here is found in the interaction between 
a forced object and a nearby wall 1 18 . The wall can be represented by an image (though not 
identical) object forced in the opposite direction 49 . As a result, the object will rotate and, 
if it is non-spherical, also glide toward the wall, as was indeed shown for a rod falling near 
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a wall 48 . Obviously, the interaction of an alignable object with a wall will also degrade the 
alignment. 

An important distinction between regular and irregular objects, which we have not dealt 
with here, concerns many-body interactions in forced systems. A pair of forced spheres 
does not develop any relative translational velocity 1 . The same holds for a pair of forced 
uniform ellipsoids to order 1 /R 3 (for an ellipsoid, the components of <f> which correspond 
to the translational velocity vanish 32 ). For a suspension of many objects this implies that 
two-body effects on relative motion are either absent (spheres) or negligible at low volume 
fraction (ellipsoids). By contrast, as we have shown here, a pair of irregular objects develops a 
relative velocity already at order 1/R 2 , which should lead to significant two-body interactions 
in a suspension. This may bring about qualitative differences between driven suspensions of 
regular and irregular objects in relation to such phenomena as sedimentation. 

This work shows that asymmetry in sedimenting objects leads to a wealth of hydrody¬ 
namic interaction effects not seen for spheres. This study was undertaken to assess how 
interactions disrupt the rotational synchronization of such objects. However it proves to 
have striking effects independent of this alignment. The prevalent repulsion, the occasional 
entrapment and the intricate quasiperiodic motions shown above are examples. These ef¬ 
fects could have significant impacts on real colloidal dispersions, e.g., in fluidized beds of 
catalyst particles. Though we have studied only pairwise interactions between identical ob¬ 
jects, many of these effects are expected to apply more generally. The general treatment of 
hydrodynamic interaction and its dependence on the shape of the interacting objects, which 
we have developed here, should prove useful in exploring these phenomena. Our work in 
progress aims to achieve a more general understanding of the rich behavior reported in Sec. 

El 

ACKNOWLEDGMENTS 

We thank Robert Deegan and Alex Leshansky for helpful discussions, and the James 
Franck Institute and Tel Aviv University for their hospitality during part of this work. This 
research has been supported by the US-Israel Binational Science Foundation (Grant no. 
2012090). 


24 


Appendix A: Notation 


The dynamics of arbitrarily shaped objects is complex and requires an elaborate notation. 
We use the following notation regarding vectors, tensors, and matrices: 

1. 3-vectors are denoted by an arrow, v, and unit 3-vectors by a hat, v. 

2 . 6-vectors are denoted by a calligraphic font, T. 

3. Matrices are marked by a blackboard-bold letter, e.g., M, where the dimension of the 
matrix is understood from the context. 

4. Tensors of rank 3 are denoted by a capital Greek letter, e.g., <f>. 

5. A set of N 3-vectors, representing N stokeslets, is denoted by a bold letter, e.g., 
v a = (vf,.. .,v%). 

6 . Subscripts with parentheses, e.g., M( 2 ), represent a term in a multipole expansion. 

7. I nxn is the n x n identity matrix. 

8 . Tensor multiplication — the dot notation,-denotes a contraction over one index. 

The double dot notation, :, denotes a contraction over two indices. Thus, given a 
tensor T of rank N and a tensor S of rank M > N, the tensors T • 5 and T : S are 
tensors of rank N + M — 2 and N + M — 4. For example, for T of rank 2 and 5 of 
rank 3, (5 • T) ikj = T is E skj and (T : E)j = T ks E skj . 

9. The matrix Y x obtained from the vector Y is defined as (Y y ') l3 = £i k jY k , such that, 
for any vector X, Y x ■ X — Y x X. 

Appendix B: Pair-mobility: Change of object origin 

Here we derive the transformation of the pair-mobility matrix under change of objects’ 
origins. Consider a new choice of origins given by R al = R a + h a and R b/ = R b + h b , and 
denote the objects’ properties with respect to the new origins with '. Following Ref. 9J the 
transformations for the generalized velocities and forces can be written as V x/ = [16x6 — 
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EX ) T ]V X an( ] jrxi _ ^ 6xg _j_ -^x^jrx f or a; = a, 6, where 


0 0 
-h ax 0 

Using [I 6x6 + B x ] _1 = [I 6x6 — B x ] we have 


I /aa M /afe 

•alba T\M/bb 


and B 6 = 


[i 6 x6 - m T ] 0 

0 [I 6x6 - (B 6 ) T ] 




[^6x6 — ®“] 0 

0 [1-6x6 — ® 6 ] 

(B.1) 


Appendix C: Properties of the tensor <h 


Below we provide a more detailed discussion regarding the tensor $ introduced in Sec. L i 
We consider its symmetries and its dependence on the choice of origin. We separate $ into 
a translational part— linear velocity response to a flow gradient, denoted by <I> tran , and a 
rotational part— angular velocity response to a flow gradient, denoted by <E> rot . We show that 
$tran is symmetric with respect to its last two indices while <f> rot has also an antisymmetric 
part which is the Levi-Civita tensor. In addition, we show that <f>tran depends on the choice 
of the object’s origin whereas <h rot does not, and derive the transformation of the former 
under change of origins. 

In order to prove the symmetry properties of we consider its transpose tensor $ r = 
which gives the force dipole around the object when subjected to external forcing, (■ rF ) = 
• T = $tran • F + <f> ro t • t . We write the force dipole as a sum of symmetric and anti¬ 
symmetric terms, | [(rP) + ( rF) T + e ■ r] = ^tran • F + <h rot • r, where e is the Levi-Civita 
tensor. The last equality implies that ( < htran)sfci is symmetric with respect to s and k and 
that the anti-symmetric part of ($ ro t ) s ki is \e s ki- 

Next we consider the transformation of <!> under change of origins. Let us assume that 
an object is given in a constant, arbitrary shear flow u{f) = S ■ r, where S is not necessarily 
a symmetric matrix. The object’s linear velocities measured about R and R! = R + h are 
V = S ■ R + $t ran : S and V' — S ■ (R + h) + $[ ran : S respectively. The tensor <I> ro t does 
not depend on the choice of origin since the angular velocity of the object is independent of 
that choice, uj = <I> ro t : S = $' ot : S. Using the relation V' = V — h x uj we find 


< an = 5 = ($tran ~ ^ ' $rot) = S ~ S ■ h. 
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In general, with analogy to Eq. (B. 1), we can write 


= [I 6x6 - (1) T ] • + A, 


(C.l) 


where 


0 0 \ I -8 is h k 11 = 1... 3 

and A iks — < 

-h x 0/ 0 , i = 4 ... 6 


Appendix D: Proofs of general properties of interaction multipoles 


Here we prove the two general results presented in Sec. |III| concerning the interaction 
multipoles. 

Multipole expansions are constructed by repeated projections (“reflections”), between 
the two objects, of the Green’s function and its derivatives. The self-blocks of the mobility 
matrix result from even projections, and the coupling blocks from odd projections. In our 
case G, the Oseen tensor, has even parity and scales as 1 /R. 

The Green’s function G itself appears only once in the expansion, in the first (1/i?) 
multipole. This is because the force monopoles acting on the particles are prescribed. This 
monopolar (odd) interaction appears only in the coupling blocks. The leading multipole 
appearing in the self-blocks is constructed by projecting the induced force dipole on object 
2 (proportional to VG) back onto object 1 (by another VG). Thus, this leading multipole 


is of 4th order, proportional to 1 /R 4 . This proves the first result in Sec. Ill Its particular 
manifestation for two spheres is well knowiP. 

Now, consider the nth multipole, proportional to 1/R n . Assume that it contains k G’s 
and n — k derivatives. Its parity is (— l) n ~ k . As explained above, for self-blocks k is even, 
and for coupling blocks it is odd. Hence, the parity of the nth multipole is (—l) n in the 
self-blocks and (—l) n+1 in the coupling blocks. This proves the second result. 


Appendix E: General Form of Mg, 

Below we provide a general form of the matrix M®^, the 2nd-order multipole of the 
coupling block in the pair-mobility matrix, and point out the number of its independent 
components. This is done by decomposing the tensors $ and 0 to their symmetric and 
anti-symmetric parts. Without loss of generality we choose the separation vector between 


27 





the two objects to be along the x axis, R = x. For two not necessarily identical objects the 
matrix has the form 
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Ad _ Ah _ Ab _ Ab 
J *~XX -^XX ^yx ^Z; 
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V 


V 


Ay x 0 0 

A-zx 0 0 / 

( T a xx 0 

T y a x 0 1 
y T? x -1 0 y 
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0 
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0 / 
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(E.l) 


where the Afj and are functions of R and the shape and orientation of object x, {x = a,b). 
For two identical (in shape and orientation) objects we have 
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